function bvp2

% code to solve the BVP
%  y'' = f(x, y, y')  for  x0 < x < x1   
%  y(x0) = y0  ,    y(x1) = y1

% get(gcf)
set(gcf,'Position', [1136 632 571 380]);

% set boundary conditions
x0 = 0.0;	y0 = 0.0;
x1 = 1.0;	y1 = 0.0;

% parameters for calculation
nx = 6;
error = 0.000001;

% start off with a linear solution
x=linspace(x0,x1,nx+2); 
gam=1;
for ix=1:nx+2
	y(ix)=-gam*x(ix)*(1-x(ix));
end;
dx=x(2)-x(1);
dxx = dx*dx;
err=1;

counter=0;
while err > error

	% calculate the coefficients of finite difference equation
	a=zeros(1,nx); c=zeros(1,nx); v=zeros(1,nx); u=zeros(1,nx);
	for j = 2:nx+1
		jj=j-1;
		z = (y(j+1) - y(j-1))/(2*dx);
		a(jj) = 2 + dxx*fy(x(j), y(j), z);
		c(jj) = -1 - 0.5*dx*fz(x(j), y(j), z);
		v(jj) = - 2*y(j) + y(j+1) + y(j-1) - dxx*f(x(j), y(j), z);
	end;

	% Newton iteration
	v(1) = v(1)/a(1);
	u(1) = - (2 + c(1))/a(1);
	for j = 2:nx
		xl = a(j) - c(j)*u(j-1);
		v(j) = (v(j) - c(j)*v(j-1))/xl;
		u(j) = - (2 + c(j))/xl;
	end;
	vv = v(nx);
	y(nx+1) = y(nx+1) + vv;
	err = abs(vv);
	for jj = nx:-1:2
		vv = v(jj-1) - u(jj-1)*vv;
		err = max(err, abs(vv));
		y(jj) = y(jj) + vv;
	end;
	counter=counter+1;
end;

newton_iterations=counter

subplot(2,1,1)
% plot computed solution
plot(x,y,'--ko','LineWidth',1,'MarkerSize',7)
hold on
set(gca,'ytick',[0  0.05 0.1 0.15]);
set(gca,'xtick',[0  0.2 0.4 0.6 0.8 1.0]);
axis([0 1 0 0.15])
xlabel('x-axis','FontSize',14,'FontWeight','bold')
ylabel('Solution','FontSize',14,'FontWeight','bold')

%plot exact solution
c0=0.7585822995;
xe=linspace(x0,x1,45);
for ix=1:45
	exact(ix)=log(2*(c0^2)*(1-tanh(c0*(xe(ix)-1/2))^2));
end;
plot(xe,exact,'k','LineWidth',1)
legend(' N = 6',' Exact',2)
% Set the fontsize to 14 for the plot
set(gca,'FontSize',14); 
% Set legend font to 14/bold                            		
set(findobj(gcf,'tag','legend'),'FontSize',14,'FontWeight','bold');  
hold off

% iteration to determine the error as number of points used is increased
ii=0;
npts=5;
nmin=log(6)/log(10);
pts=linspace(nmin,2.98,npts);
for i=1:npts
	
	ii=ii+1; points(ii)=round(10^pts(i))
	nx=points(ii);
	error = 0.0000000001;

	% start off with a linear solution
	x=linspace(x0,x1,nx+2);
	for ix=1:nx+2
		y(ix)=-gam*x(ix)*(1-x(ix));
		exact2(ix)=log(2*(c0^2)*(1-tanh(c0*(x(ix)-1/2))^2));
	end;
	dx=x(2)-x(1);
	dxx = dx*dx;
	err=1;
	while err > error

	% calculate the coefficients of finite difference equation
	a=zeros(1,nx); c=zeros(1,nx); v=zeros(1,nx); u=zeros(1,nx);
	for j = 2:nx+1
		jj=j-1;
		z = (y(j+1) - y(j-1))/(2*dx);
		a(jj) = 2 + dxx*fy(x(j), y(j), z);
		c(jj) = -1 - 0.5*dx*fz(x(j), y(j), z);
		v(jj) = - 2*y(j) + y(j+1) + y(j-1) - dxx*f(x(j), y(j), z);
	end;

	% Newton iteration
	v(1) = v(1)/a(1);
	u(1) = - (2 + c(1))/a(1);
	for j = 2:nx
		xl = a(j) - c(j)*u(j-1);
		v(j) = (v(j) - c(j)*v(j-1))/xl;
		u(j) = - (2 + c(j))/xl;
	end;
	vv = v(nx);
	y(nx+1) = y(nx+1) + vv;
	err = abs(vv);
	for jj = nx:-1:2
		vv = v(jj-1) - u(jj-1)*vv;
		err = max(err, abs(vv));
		y(jj) = y(jj) + vv;
	end;

	end;
	errorM(ii)=norm(exact2-y,inf);
end;

subplot(2,1,2)
loglog(points,errorM,'--ok','LineWidth',1,'MarkerSize',7)
hold on
set(gca,'ytick',[1e-9 1e-7 1e-5 1e-3]);
axis([1 1e3 1e-9 1e-3])
grid on
set(gca,'MinorGridLineStyle','none')
xlabel('Number of Grid Points','FontSize',14,'FontWeight','bold')
ylabel('Error','FontSize',14,'FontWeight','bold')
% Set the fontsize to 14 for the plot
set(gca,'FontSize',14);
hold off


function g=f(x,y,z)
g=-exp(y);

function g=fy(x,y,z)
g=-exp(y);

function g=fz(x,y,z)
g=0;







